Testing hypotheses for pairs of species using ENMTools
Building species
Just in case we don’t have our data saved from the last round:
library(rgbif)
library(ENMTools)
ibm <- occ_search(scientificName = "Iberolacerta monticola", limit = 1500)
ibc <- occ_search(scientificName = "Iberolacerta cyreni", limit = 1500)
keep.cols <- c("species", "decimalLatitude", "decimalLongitude")
ibm <- ibm$data[,keep.cols]
ibm <- as.data.frame(unique(ibm))
ibm <- ibm[complete.cases(ibm),]
colnames(ibm) <- c("species", "lat", "lon")
ibc <- ibc$data[,keep.cols]
ibc <- as.data.frame(unique(ibc))
ibc <- ibc[complete.cases(ibc),]
colnames(ibc) <- c("species", "lat", "lon")
all.worldclim <- raster::getData(name = "worldclim", download = TRUE, res = 10, var = "bio")
spain.worldclim <- crop(all.worldclim, extent(-10, 4, 35, 45))
plot(spain.worldclim[["bio1"]])
points(ibc[,c("lon", "lat")], pch = 16, cex = 0.5, col = "red")
points(ibm[,c("lon", "lat")], pch = 16, cex = 0.5, col = "blue")
Okay, let’s turn both of these into enmtools.species objects
monticola <- enmtools.species()
monticola$presence.points <- ibm[,c("lon", "lat")]
monticola$species.name <- "Iberolacerta monticola"
monticola$range <- background.raster.buffer(monticola$presence.points, 50000, mask = spain.worldclim)
plot(monticola)
cyreni <- enmtools.species()
cyreni$presence.points <- ibc[,c("lon", "lat")]
cyreni$species.name <- "Iberolacerta cyreni"
cyreni$range <- background.raster.buffer(cyreni$presence.points, 50000, mask = spain.worldclim)
plot(cyreni)
cyreni <- check.species(cyreni)
monticola <- check.species(monticola)Now let’s build a quick ENM for each species. Since the background for cyreni is so small, we’ve got to restrict ourselves to only using 100 background points.
monticola.glm <- enmtools.glm(monticola, spain.worldclim, f = pres ~ poly(bio1, 4) + poly(bio8, 4), nback = 100, test.prop = 0.3)
cyreni.glm <- enmtools.glm(cyreni, spain.worldclim, f = pres ~ poly(bio1, 4) + poly(bio8, 4), nback = 100, test.prop = 0.3)
cyreni.glm##
##
## Formula: presence ~ poly(bio1, 4) + poly(bio8, 4)
## <environment: 0x7f91870674d0>
##
##
## Data table (top ten lines):
##
## | | Longitude| Latitude| bio1| bio2| bio3| bio4| bio5| bio6| bio7| bio8| bio9| bio10| bio11| bio12| bio13| bio14| bio15| bio16| bio17| bio18| bio19| presence|
## |:--|---------:|--------:|----:|----:|----:|----:|----:|----:|----:|----:|----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|--------:|
## |1 | -5.298285| 40.30419| 89| 111| 37| 6278| 267| -29| 296| 22| 174| 175| 16| 646| 74| 19| 34| 206| 80| 84| 187| 1|
## |2 | -5.268747| 40.20154| 89| 111| 37| 6278| 267| -29| 296| 22| 174| 175| 16| 646| 74| 19| 34| 206| 80| 84| 187| 1|
## |3 | -5.045487| 40.23004| 115| 113| 37| 6357| 293| -8| 301| 99| 200| 202| 40| 483| 58| 13| 34| 153| 61| 62| 133| 1|
## |6 | -3.818184| 40.99835| 94| 105| 37| 6255| 262| -21| 283| 76| 179| 179| 22| 558| 65| 20| 28| 169| 84| 92| 144| 1|
## |7 | -5.276188| 40.23706| 89| 111| 37| 6278| 267| -29| 296| 22| 174| 175| 16| 646| 74| 19| 34| 206| 80| 84| 187| 1|
## |8 | -3.826836| 40.85220| 94| 105| 37| 6255| 262| -21| 283| 76| 179| 179| 22| 558| 65| 20| 28| 169| 84| 92| 144| 1|
## |9 | -3.817753| 40.94713| 94| 105| 37| 6255| 262| -21| 283| 76| 179| 179| 22| 558| 65| 20| 28| 169| 84| 92| 144| 1|
## |11 | -3.869114| 40.98392| 80| 106| 37| 6102| 249| -32| 281| 94| 164| 164| 11| 618| 74| 24| 27| 190| 96| 106| 157| 1|
## |12 | -5.232444| 40.27478| 89| 111| 37| 6278| 267| -29| 296| 22| 174| 175| 16| 646| 74| 19| 34| 206| 80| 84| 187| 1|
## |13 | -5.023499| 40.44747| 85| 111| 37| 6218| 262| -33| 295| 46| 169| 170| 13| 605| 74| 19| 32| 189| 81| 86| 164| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2490 -0.4420 -0.1205 0.3031 2.0966
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6183 1.1693 0.529 0.59694
## poly(bio1, 4)1 -45.3305 19.4085 -2.336 0.01951 *
## poly(bio1, 4)2 -21.5334 20.2501 -1.063 0.28761
## poly(bio1, 4)3 -21.4025 13.2599 -1.614 0.10651
## poly(bio1, 4)4 -15.1579 9.6291 -1.574 0.11545
## poly(bio8, 4)1 -49.5074 22.2360 -2.226 0.02598 *
## poly(bio8, 4)2 41.5064 24.8239 1.672 0.09452 .
## poly(bio8, 4)3 -40.4693 14.3816 -2.814 0.00489 **
## poly(bio8, 4)4 13.7627 10.9005 1.263 0.20674
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 185.763 on 166 degrees of freedom
## Residual deviance: 85.847 on 158 degrees of freedom
## AIC: 125.33
##
## Number of Fisher Scoring iterations: 9
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 67
## n absences : 100
## AUC : 0.9367164
## cor : 0.5772312
## max TPR+TNR at : -0.1505364
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 67
## n absences : 10000
## AUC : 0.7931597
## cor : 0.08500655
## max TPR+TNR at : 0.4623617
##
##
## Proportion of data wittheld for model fitting: 0.3
##
## Model fit (test data): class : ModelEvaluation
## n presences : 29
## n absences : 100
## AUC : 0.9136207
## cor : 0.5852398
## max TPR+TNR at : 0.7135983
##
##
## Environment space model fit (test data): class : ModelEvaluation
## n presences : 29
## n absences : 10000
## AUC : 0.7843897
## cor : 0.05282975
## max TPR+TNR at : 0.6711178
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 2.220446e-16, 1 (min, max)
##
##
##
## Notes:
monticola.glm##
##
## Formula: presence ~ poly(bio1, 4) + poly(bio8, 4)
## <environment: 0x7f9188351e10>
##
##
## Data table (top ten lines):
##
## | | Longitude| Latitude| bio1| bio2| bio3| bio4| bio5| bio6| bio7| bio8| bio9| bio10| bio11| bio12| bio13| bio14| bio15| bio16| bio17| bio18| bio19| presence|
## |:--|---------:|--------:|----:|----:|----:|----:|----:|----:|----:|----:|----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|--------:|
## |2 | -7.188726| 43.43051| 128| 84| 41| 4080| 239| 38| 201| 83| 180| 182| 78| 959| 124| 33| 34| 345| 126| 142| 310| 1|
## |6 | -6.728039| 43.55888| 137| 83| 42| 3957| 244| 48| 196| 115| 187| 190| 89| 852| 114| 33| 31| 298| 124| 139| 253| 1|
## |7 | -6.823801| 43.39811| 121| 88| 41| 4237| 239| 28| 211| 97| 175| 178| 70| 920| 121| 36| 31| 320| 135| 148| 279| 1|
## |8 | -6.732560| 43.54083| 137| 83| 42| 3957| 244| 48| 196| 115| 187| 190| 89| 852| 114| 33| 31| 298| 124| 139| 253| 1|
## |9 | -6.650406| 43.40437| 117| 90| 42| 4292| 237| 23| 214| 93| 172| 174| 65| 926| 123| 38| 30| 320| 140| 153| 274| 1|
## |11 | -6.910225| 43.29794| 117| 91| 41| 4472| 240| 20| 220| 67| 175| 176| 63| 942| 122| 34| 33| 334| 129| 142| 300| 1|
## |13 | -6.817537| 43.28005| 105| 93| 41| 4614| 233| 7| 226| 54| 165| 166| 49| 983| 127| 38| 32| 345| 141| 152| 309| 1|
## |14 | -6.831939| 43.38955| 121| 88| 41| 4237| 239| 28| 211| 97| 175| 178| 70| 920| 121| 36| 31| 320| 135| 148| 279| 1|
## |16 | -6.499937| 43.47067| 117| 91| 42| 4248| 237| 23| 214| 94| 171| 174| 66| 913| 122| 40| 29| 313| 145| 156| 261| 1|
## |17 | -6.619466| 43.42050| 117| 90| 42| 4292| 237| 23| 214| 93| 172| 174| 65| 926| 123| 38| 30| 320| 140| 153| 274| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2300 -0.4855 0.9284 1.1686 3.3594
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.08128 0.13634 -0.596 0.551097
## poly(bio1, 4)1 -3.31977 4.60674 -0.721 0.471136
## poly(bio1, 4)2 -36.10257 12.21081 -2.957 0.003110 **
## poly(bio1, 4)3 14.07981 10.20146 1.380 0.167533
## poly(bio1, 4)4 -26.55105 7.86313 -3.377 0.000734 ***
## poly(bio8, 4)1 -5.56696 2.46036 -2.263 0.023656 *
## poly(bio8, 4)2 1.69600 2.38727 0.710 0.477436
## poly(bio8, 4)3 -5.18533 2.49142 -2.081 0.037409 *
## poly(bio8, 4)4 -7.14763 2.59940 -2.750 0.005964 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 695.92 on 350 degrees of freedom
## Residual deviance: 593.46 on 342 degrees of freedom
## AIC: 670.46
##
## Number of Fisher Scoring iterations: 8
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 251
## n absences : 100
## AUC : 0.7297211
## cor : 0.2972758
## max TPR+TNR at : -0.07814603
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 251
## n absences : 10000
## AUC : 0.7010849
## cor : 0.119513
## max TPR+TNR at : 0.4534149
##
##
## Proportion of data wittheld for model fitting: 0.3
##
## Model fit (test data): class : ModelEvaluation
## n presences : 107
## n absences : 100
## AUC : 0.7116355
## cor : 0.2518333
## max TPR+TNR at : -0.058809
##
##
## Environment space model fit (test data): class : ModelEvaluation
## n presences : 107
## n absences : 10000
## AUC : 0.6679748
## cor : 0.06633703
## max TPR+TNR at : 0.2257457
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 2.220446e-16, 0.8932101 (min, max)
##
##
##
## Notes:
Breadth and overlap
Okay, now let’s look at breadth of species ENMs in geographic space.
raster.breadth(cyreni.glm)## $B1
## [1] 0.8788605
##
## $B2
## [1] 0.3010483raster.breadth(monticola.glm)## $B1
## [1] 0.9059536
##
## $B2
## [1] 0.4387881How about in environment space?
env.breadth(cyreni.glm, spain.worldclim)## $env.B2
## [1] 0.3984402
##
## $B2.plot
env.breadth(monticola.glm, spain.worldclim)## $env.B2
## [1] 0.3340117
##
## $B2.plot
Now let’s calculate niche overlap in geographic space and environment space.
raster.overlap(monticola.glm, cyreni.glm)## $D
## [1] 0.6977142
##
## $I
## [1] 0.8986253
##
## $rank.cor
## [1] 0.8570302env.overlap(monticola.glm, cyreni.glm, spain.worldclim)## $env.D
## [1] 0.5388819
##
## $env.I
## [1] 0.6880965
##
## $env.cor
## [1] 0.4124475
##
## $env.D.plot
##
## $env.I.plot
##
## $env.cor.plot
So the models are much more similar in geographic space than they are in environment space. What could this mean? DISCUSS
Identity test
Let’s do an identity/equivalency test!
id.glm <- identity.test(monticola, cyreni, spain.worldclim, type = "glm", f = pres ~ poly(bio1, 4) + poly(bio8, 4), nreps = 20, nback = 100)
id.glm##
##
##
##
## Identity test Iberolacerta monticola vs. Iberolacerta cyreni
##
## objectentity test p-values:
## D I rank.cor env.D env.I env.cor
## 0.04761905 0.09523810 0.23809524 0.23809524 0.23809524 0.23809524
##
##
## Replicates:
##
##
## | | D| I| rank.cor| env.D| env.I| env.cor|
## |:---------|---------:|---------:|---------:|---------:|---------:|---------:|
## |empirical | 0.7518235| 0.9469363| 0.9295694| 0.7298343| 0.8900250| 0.8846858|
## |rep 1 | 0.9026713| 0.9872296| 0.9853723| 0.9183663| 0.9807259| 0.9605894|
## |rep 2 | 0.8831265| 0.9743586| 0.9485421| 0.8702099| 0.9628094| 0.9213665|
## |rep 3 | 0.8768048| 0.9761098| 0.9555575| 0.8959364| 0.9753094| 0.9483066|
## |rep 4 | 0.9212783| 0.9898480| 0.9893446| 0.9057748| 0.9767442| 0.9683371|
## |rep 5 | 0.9021508| 0.9823899| 0.9647467| 0.9155530| 0.9807764| 0.9616315|
Background test
Let’s do a background test!
bg.glm <- background.test(monticola, cyreni, spain.worldclim, type = "glm", f = pres ~ poly(bio1, 4) + poly(bio8, 4), nreps = 20, nback = 100, test.type = "symmetric")## D I rank.cor env.D env.I env.cor
## empirical 0.7847587 0.9602411 0.95839412 0.7420224 0.8957305 0.900453878
## rep 1 0.8651889 0.9838006 -0.34402277 0.7436440 0.9270632 -0.331782392
## rep 2 0.7865694 0.9166384 -0.11171020 0.4324687 0.6928510 -0.436137060
## rep 3 0.9097958 0.9888669 0.65496021 0.8025011 0.9415045 0.609072867
## rep 4 0.8820885 0.9827633 0.78527174 0.8021663 0.9580440 0.761307159
## rep 5 0.8799236 0.9854709 0.27812430 0.7682087 0.9185553 0.462282068
## rep 6 0.9322644 0.9871165 0.78753599 0.6418691 0.8620046 0.227469968
## rep 7 0.6147337 0.8369891 -0.54974959 0.5627812 0.7547207 -0.117212142
## rep 8 0.8339719 0.9688869 -0.79254629 0.6132130 0.8331794 -0.420741925
## rep 9 0.9009113 0.9869945 0.21905689 0.8947915 0.9778109 0.437667107
## rep 10 0.8599917 0.9774707 -0.16457533 0.6313168 0.8094308 -0.005796276
## rep 11 0.8429407 0.9770776 0.02337592 0.6287944 0.8638191 -0.391687540
## rep 12 0.8936726 0.9856197 0.14287734 0.6955805 0.9086622 -0.371173593
## rep 13 0.7011788 0.8971583 -0.27785201 0.7709377 0.9194123 0.181089319
## rep 14 0.7498888 0.8855333 -0.26720290 0.4683470 0.6970489 -0.521037512
## rep 15 0.8870873 0.9759267 0.60554310 0.7477945 0.9077167 0.043933447
## rep 16 0.8834884 0.9842235 0.66730122 0.7687547 0.9240908 0.418878461
## rep 17 0.8969396 0.9870014 0.07937109 0.7181077 0.9103742 -0.596978113
## rep 18 0.8105311 0.9517275 -0.05681624 0.8117832 0.9461840 0.441746975
## rep 19 0.9120719 0.9929056 0.36416909 0.8532650 0.9802262 0.630160268
## rep 20 0.9102264 0.9857067 0.61854130 0.7072515 0.8990217 0.523940015bg.glm## D I rank.cor env.D env.I env.cor
## 0.3809524 0.5714286 0.0952381 1.0476190 0.7619048 0.0952381
##
##
## | | D| I| rank.cor| env.D| env.I| env.cor|
## |:---------|---------:|---------:|----------:|---------:|---------:|----------:|
## |empirical | 0.7847587| 0.9602411| 0.9583941| 0.7420224| 0.8957305| 0.9004539|
## |rep 1 | 0.8651889| 0.9838006| -0.3440228| 0.7436440| 0.9270632| -0.3317824|
## |rep 2 | 0.7865694| 0.9166384| -0.1117102| 0.4324687| 0.6928510| -0.4361371|
## |rep 3 | 0.9097958| 0.9888669| 0.6549602| 0.8025011| 0.9415045| 0.6090729|
## |rep 4 | 0.8820885| 0.9827633| 0.7852717| 0.8021663| 0.9580440| 0.7613072|
## |rep 5 | 0.8799236| 0.9854709| 0.2781243| 0.7682087| 0.9185553| 0.4622821|
Ecospat
ib.ecospat.id <- enmtools.ecospat.id(monticola, cyreni, spain.worldclim, nreps = 100, layers = c("bio1", "bio8"))
ib.ecospat.id##
##
##
##
## Ecospat identity test Iberolacerta monticola vs. Iberolacerta cyreni
##
## ecospat.id test empirical overlaps:
## $D
## [1] 0.05644239
##
## $I
## [1] 0.2219478
##
##
##
## ecospat.id test p-values:
## D I
## 0 0
## NULLib.ecospat.bg <- enmtools.ecospat.bg(monticola, cyreni, spain.worldclim, nreps = 100, layers = c("bio1", "bio8"))
ib.ecospat.bg##
##
##
##
## Ecospat background test asymmetric Iberolacerta monticola vs. Iberolacerta cyreni
##
## ecospat.bg test empirical overlaps:
## $D
## [1] 0.05644239
##
## $I
## [1] 0.2219478
##
##
##
## ecospat.bg test p-values:
## D I
## 0.01980198 0.01980198
## NULLMOSES
ib.moses <- moses.list(list(cyreni, monticola), spain.worldclim[[c(1,8)]])
ib.moses## $separate.glms
## $separate.glms$`Iberolacerta cyreni`
##
##
## Formula: presence ~ bio1 + bio8
## <environment: 0x7f917db3d3f0>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio8| presence|
## |---------:|--------:|----:|----:|--------:|
## | -5.298285| 40.30419| 89| 22| 1|
## | -5.268747| 40.20154| 89| 22| 1|
## | -5.045487| 40.23004| 115| 99| 1|
## | -5.373484| 40.27068| 93| 26| 1|
## | -5.329139| 40.22420| 89| 22| 1|
## | -3.818184| 40.99835| 94| 76| 1|
## | -5.276188| 40.23706| 89| 22| 1|
## | -3.826836| 40.85220| 94| 76| 1|
## | -3.817753| 40.94713| 94| 76| 1|
## | -3.967111| 40.87878| 80| 94| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6632 -0.6795 -0.2751 0.5057 2.3797
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.471994 1.350645 7.013 2.33e-12 ***
## bio1 -0.081621 0.012490 -6.535 6.37e-11 ***
## bio8 -0.009363 0.005782 -1.619 0.105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 266.17 on 249 degrees of freedom
## Residual deviance: 182.45 on 247 degrees of freedom
## AIC: 241.76
##
## Number of Fisher Scoring iterations: 5
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 96
## n absences : 154
## AUC : 0.8552151
## cor : 0.5809062
## max TPR+TNR at : 0.7695674
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 96
## n absences : 10000
## AUC : 0.8121542
## cor : 0.1117736
## max TPR+TNR at : 0.7410556
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0.0006689058, 0.99925 (min, max)
##
##
##
## Notes:
##
## $separate.glms$`Iberolacerta monticola`
##
##
## Formula: presence ~ bio1 + bio8
## <environment: 0x7f918f051230>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio8| presence|
## |---------:|--------:|----:|----:|--------:|
## | -6.699292| 43.45172| 121| 97| 1|
## | -7.188726| 43.43051| 128| 83| 1|
## | -6.703289| 43.46805| 121| 97| 1|
## | -6.504768| 43.51037| 134| 113| 1|
## | -6.556106| 43.52801| 134| 113| 1|
## | -6.728039| 43.55888| 137| 115| 1|
## | -6.823801| 43.39811| 121| 97| 1|
## | -6.732560| 43.54083| 137| 115| 1|
## | -6.650406| 43.40437| 117| 93| 1|
## | -6.541960| 43.45676| 117| 93| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3221 -0.9727 -0.8129 1.0716 1.5583
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.442609 0.384387 3.753 0.000175 ***
## bio1 -0.005827 0.004177 -1.395 0.162945
## bio8 -0.010681 0.003486 -3.064 0.002183 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 992.59 on 877 degrees of freedom
## Residual deviance: 966.99 on 875 degrees of freedom
## AIC: 1191.7
##
## Number of Fisher Scoring iterations: 4
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 358
## n absences : 520
## AUC : 0.6045713
## cor : 0.1852136
## max TPR+TNR at : 0.2538606
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 358
## n absences : 10000
## AUC : 0.5311218
## cor : 0.01901879
## max TPR+TNR at : 0.4248563
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0.1508468, 0.7595605 (min, max)
##
##
##
## Notes:
##
##
## $combined.glm
##
##
## Formula: presence ~ bio1 + bio8
## <environment: 0x7f9178ea6240>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio8| presence|
## |---------:|--------:|----:|----:|--------:|
## | -5.298285| 40.30419| 89| 22| 1|
## | -5.268747| 40.20154| 89| 22| 1|
## | -5.045487| 40.23004| 115| 99| 1|
## | -5.373484| 40.27068| 93| 26| 1|
## | -5.329139| 40.22420| 89| 22| 1|
## | -3.818184| 40.99835| 94| 76| 1|
## | -5.276188| 40.23706| 89| 22| 1|
## | -3.826836| 40.85220| 94| 76| 1|
## | -3.817753| 40.94713| 94| 76| 1|
## | -3.967111| 40.87878| 80| 94| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6261 -0.9435 -0.7464 1.0307 1.5630
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.505939 0.364621 6.873 6.30e-12 ***
## bio1 -0.016862 0.003726 -4.525 6.03e-06 ***
## bio8 -0.008527 0.002816 -3.028 0.00246 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1258.8 on 1127 degrees of freedom
## Residual deviance: 1196.6 on 1125 degrees of freedom
## AIC: 1493.2
##
## Number of Fisher Scoring iterations: 4
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 454
## n absences : 674
## AUC : 0.6510592
## cor : 0.2531627
## max TPR+TNR at : 0.2726938
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 454
## n absences : 10000
## AUC : 0.4603159
## cor : -0.02365221
## max TPR+TNR at : 0.3040349
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 60, 84, 5040 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 4, 35, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0.09451035, 0.8595204 (min, max)
##
##
##
## Notes:
##
## $separate.aic
## [1] 1433.427
##
## $combined.aic
## [1] 1493.171Okay not much exciting there. Let’s load in an enmtools.clade object though.
data(iberolacerta.clade)
iberolacerta.clade##
##
## An enmtools.clade object with 7 species
##
## Species names:
## monticola martinezricai cyreni horvathi aurelioi aranica bonnali
##
## Tree:
##
## Phylogenetic tree with 7 tips and 6 internal nodes.
##
## Tip labels:
## monticola, martinezricai, cyreni, horvathi, aurelioi, aranica, ...
##
## Rooted; includes branch lengths.
##
##
## Data Summary:
##
##
## | |species.names |in.tree |presence |background |range |
## |:-------------|:-------------|:-------|:--------|:----------|:-------|
## |monticola |monticola |TRUE |260 |0 |present |
## |martinezricai |martinezricai |TRUE |3 |0 |present |
## |cyreni |cyreni |TRUE |76 |0 |present |
## |horvathi |horvathi |TRUE |4 |0 |present |
## |aurelioi |aurelioi |TRUE |28 |0 |present |
## |aranica |aranica |TRUE |14 |0 |present |
## |bonnali |bonnali |TRUE |56 |0 |present |plot(iberolacerta.clade$tree)
all.worldclim <- raster::getData(name = "worldclim", download = TRUE, res = 10, var = "bio")
euro.worldclim <- crop(all.worldclim, extent(-10, 17, 39, 48))
euro.worldclim <- euro.worldclim[[c(1,7,12,14)]]
raster.cor.matrix(env = euro.worldclim)## bio1 bio7 bio12 bio14
## bio1 1.00000000 -0.05329392 -0.6084667 -0.72871981
## bio7 -0.05329392 1.00000000 -0.2913319 -0.07381576
## bio12 -0.60846674 -0.29133189 1.0000000 0.77619752
## bio14 -0.72871981 -0.07381576 0.7761975 1.00000000ib.moses.1 <- moses.list(list(iberolacerta.clade$species$aurelioi, iberolacerta.clade$species$aranica), env = euro.worldclim)
ib.moses.1## $separate.glms
## $separate.glms$aurelioi
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f917719e080>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | 1.568714| 42.82648| 69| 235| 1087| 74| 1|
## | 1.473000| 42.50300| 31| 225| 1322| 88| 1|
## | 1.492613| 42.77308| 53| 232| 1191| 82| 1|
## | -1.480000| 42.58000| 113| 265| 850| 42| 1|
## | 1.230000| 42.75000| 58| 234| 1170| 80| 1|
## | 1.230000| 42.66000| 60| 234| 1142| 77| 1|
## | 1.350000| 42.67000| 53| 232| 1191| 82| 1|
## | -1.470000| 42.67000| 111| 260| 977| 48| 1|
## | 1.360000| 42.58000| 31| 225| 1322| 88| 1|
## | 1.291700| 42.61700| 60| 234| 1142| 77| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9669 -0.5374 -0.3234 -0.1410 2.5799
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.636149 13.759368 -0.555 0.579
## bio1 -0.018214 0.043207 -0.422 0.673
## bio7 0.018223 0.044414 0.410 0.682
## bio12 0.005748 0.005520 1.041 0.298
## bio14 -0.019642 0.091601 -0.214 0.830
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 77.632 on 156 degrees of freedom
## Residual deviance: 58.807 on 152 degrees of freedom
## AIC: 39.535
##
## Number of Fisher Scoring iterations: 5
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 28
## n absences : 129
## AUC : 0.8322259
## cor : 0.4120869
## max TPR+TNR at : 0.5867616
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 28
## n absences : 10000
## AUC : 0.7395179
## cor : 0.04517637
## max TPR+TNR at : 0.6425447
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0.004775502, 0.9965576 (min, max)
##
##
##
## Notes:
##
## $separate.glms$aranica
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f91783e78c8>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | 1.572401| 42.99966| 105| 249| 871| 57| 1|
## | 1.103258| 42.77680| 41| 229| 1271| 86| 1|
## | 0.860000| 42.75000| 37| 228| 1294| 87| 1|
## | 0.860000| 42.84000| 82| 244| 1051| 72| 1|
## | -0.980000| 42.84000| 74| 243| 1096| 63| 1|
## | 0.980000| 42.75000| 37| 228| 1294| 87| 1|
## | 0.923200| 42.70100| 37| 228| 1294| 87| 1|
## | 0.804300| 42.60900| 33| 228| 1294| 86| 1|
## | 0.801100| 42.69900| 59| 237| 1166| 78| 1|
## | 0.798000| 42.78900| 59| 237| 1166| 78| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.73516 -0.34933 -0.21089 -0.09113 1.93524
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -92.259128 68.309702 -1.351 0.1768
## bio1 0.049301 0.076988 0.640 0.5219
## bio7 0.250015 0.226881 1.102 0.2705
## bio12 0.005322 0.011581 0.460 0.6459
## bio14 0.324348 0.178154 1.821 0.0687 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 38.816 on 117 degrees of freedom
## Residual deviance: 24.374 on 113 degrees of freedom
## AIC: 22.193
##
## Number of Fisher Scoring iterations: 5
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 14
## n absences : 104
## AUC : 0.8839286
## cor : 0.4209443
## max TPR+TNR at : 1.151452
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 14
## n absences : 10000
## AUC : 0.65425
## cor : 0.02759504
## max TPR+TNR at : 0.7596942
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 2.220446e-16, 1 (min, max)
##
##
##
## Notes:
##
##
## $combined.glm
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9175236ca0>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | 1.568714| 42.82648| 69| 235| 1087| 74| 1|
## | 1.473000| 42.50300| 31| 225| 1322| 88| 1|
## | 1.492613| 42.77308| 53| 232| 1191| 82| 1|
## | -1.480000| 42.58000| 113| 265| 850| 42| 1|
## | 1.230000| 42.75000| 58| 234| 1170| 80| 1|
## | 1.230000| 42.66000| 60| 234| 1142| 77| 1|
## | 1.350000| 42.67000| 53| 232| 1191| 82| 1|
## | -1.470000| 42.67000| 111| 260| 977| 48| 1|
## | 1.360000| 42.58000| 31| 225| 1322| 88| 1|
## | 1.291700| 42.61700| 60| 234| 1142| 77| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8222 -0.5040 -0.3103 -0.1606 2.8669
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.187e+01 1.521e+01 -0.780 0.435
## bio1 -3.498e-04 3.187e-02 -0.011 0.991
## bio7 2.175e-02 4.950e-02 0.439 0.660
## bio12 3.568e-03 4.422e-03 0.807 0.420
## bio14 4.277e-02 6.923e-02 0.618 0.537
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 116.449 on 274 degrees of freedom
## Residual deviance: 90.152 on 270 degrees of freedom
## AIC: 56.04
##
## Number of Fisher Scoring iterations: 5
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 42
## n absences : 233
## AUC : 0.8236767
## cor : 0.3666748
## max TPR+TNR at : 0.5650646
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 42
## n absences : 10000
## AUC : 0.8434238
## cor : 0.09020671
## max TPR+TNR at : 0.6375467
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0.001156367, 0.9986404 (min, max)
##
##
##
## Notes:
##
## $separate.aic
## [1] 61.72772
##
## $combined.aic
## [1] 56.03959ib.moses.2 <- moses.list(list(iberolacerta.clade$species$monticola, iberolacerta.clade$species$martinezricai), env = euro.worldclim)
ib.moses.2## $separate.glms
## $separate.glms$monticola
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f91672cb8e8>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -5.171215| 43.06957| 78| 249| 917| 48| 1|
## | -6.036635| 43.02531| 76| 246| 1012| 46| 1|
## | -7.679727| 40.38852| 137| 247| 1143| 11| 1|
## | -7.790437| 40.30959| 129| 242| 1231| 12| 1|
## | -7.473340| 43.78935| 140| 179| 931| 32| 1|
## | -6.575039| 42.91070| 84| 247| 1012| 39| 1|
## | -5.132756| 43.49572| 133| 190| 822| 40| 1|
## | -7.787378| 40.39362| 137| 247| 1143| 11| 1|
## | -4.941888| 43.35310| 128| 194| 843| 41| 1|
## | -7.621731| 40.34170| 101| 229| 1514| 15| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6336 -0.8287 -0.5222 0.8565 2.4616
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 24.0219102 4.3648202 5.504 3.72e-08 ***
## bio1 -0.0782058 0.0125768 -6.218 5.03e-10 ***
## bio7 -0.0498766 0.0087542 -5.697 1.22e-08 ***
## bio12 -0.0013751 0.0007543 -1.823 0.0683 .
## bio14 -0.0755313 0.0186510 -4.050 5.13e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 720.87 on 771 degrees of freedom
## Residual deviance: 621.74 on 767 degrees of freedom
## AIC: 945.6
##
## Number of Fisher Scoring iterations: 4
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 260
## n absences : 512
## AUC : 0.7444561
## cor : 0.3552669
## max TPR+TNR at : -0.01515553
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 260
## n absences : 10000
## AUC : 0.4879765
## cor : 0.01871728
## max TPR+TNR at : 0.3657356
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0.0004582117, 0.9997647 (min, max)
##
##
##
## Notes:
##
## $separate.glms$martinezricai
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f918f98f3c8>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -6.047705| 40.52360| 113| 287| 610| 11| 1|
## | -6.230000| 40.48000| 124| 284| 635| 10| 1|
## | -6.110000| 40.48000| 134| 291| 529| 9| 1|
## | -6.250000| 40.91667| 122| 286| 569| 10| 0|
## | -6.083333| 40.91667| 119| 289| 518| 10| 0|
## | -5.916667| 40.91667| 120| 293| 463| 10| 0|
## | -6.583333| 40.75000| 124| 275| 716| 10| 0|
## | -6.416667| 40.75000| 123| 280| 649| 10| 0|
## | -6.250000| 40.75000| 122| 284| 587| 10| 0|
## | -6.083333| 40.75000| 120| 289| 541| 10| 0|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4859 -0.3534 -0.3025 -0.2755 1.2744
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -62.57606 147.32448 -0.425 0.671
## bio1 -0.09859 0.28262 -0.349 0.727
## bio7 0.25413 0.50224 0.506 0.613
## bio12 0.01804 0.03804 0.474 0.635
## bio14 -0.88155 2.14362 -0.411 0.681
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 8.3178 on 39 degrees of freedom
## Residual deviance: 8.0123 on 35 degrees of freedom
## AIC: 13.926
##
## Number of Fisher Scoring iterations: 4
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 3
## n absences : 37
## AUC : 0.6711712
## cor : 0.09935648
## max TPR+TNR at : -0.2251678
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 3
## n absences : 10000
## AUC : 0.6122333
## cor : 0.005300873
## max TPR+TNR at : 0.4438694
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 2.220446e-16, 0.9999929 (min, max)
##
##
##
## Notes:
##
##
## $combined.glm
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f91773fca40>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -5.171215| 43.06957| 78| 249| 917| 48| 1|
## | -6.036635| 43.02531| 76| 246| 1012| 46| 1|
## | -7.679727| 40.38852| 137| 247| 1143| 11| 1|
## | -7.790437| 40.30959| 129| 242| 1231| 12| 1|
## | -7.473340| 43.78935| 140| 179| 931| 32| 1|
## | -6.575039| 42.91070| 84| 247| 1012| 39| 1|
## | -5.132756| 43.49572| 133| 190| 822| 40| 1|
## | -7.787378| 40.39362| 137| 247| 1143| 11| 1|
## | -4.941888| 43.35310| 128| 194| 843| 41| 1|
## | -7.621731| 40.34170| 101| 229| 1514| 15| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6204 -0.8071 -0.5020 0.8338 2.4545
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 23.6059952 4.3877943 5.380 7.45e-08 ***
## bio1 -0.0767323 0.0125613 -6.109 1.00e-09 ***
## bio7 -0.0492759 0.0087995 -5.600 2.15e-08 ***
## bio12 -0.0013220 0.0007644 -1.729 0.083747 .
## bio14 -0.0720917 0.0187971 -3.835 0.000125 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 729.19 on 811 degrees of freedom
## Residual deviance: 622.26 on 807 degrees of freedom
## AIC: 308.5
##
## Number of Fisher Scoring iterations: 5
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 263
## n absences : 549
## AUC : 0.7529106
## cor : 0.3642898
## max TPR+TNR at : 0.0441889
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 263
## n absences : 10000
## AUC : 0.4835658
## cor : 0.01907014
## max TPR+TNR at : 0.3718913
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0.000590172, 0.9997672 (min, max)
##
##
##
## Notes:
##
## $separate.aic
## [1] 959.5283
##
## $combined.aic
## [1] 308.5021montmart <- combine.species(list(iberolacerta.clade$species$monticola, iberolacerta.clade$species$martinezricai))
ib.moses.3 <- moses.list(list(montmart, iberolacerta.clade$species$cyreni), env = euro.worldclim)
ib.moses.3## $separate.glms
## $separate.glms$`monticola + martinezricai`
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9167633110>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -5.171215| 43.06957| 78| 249| 917| 48| 1|
## | -6.036635| 43.02531| 76| 246| 1012| 46| 1|
## | -7.679727| 40.38852| 137| 247| 1143| 11| 1|
## | -7.790437| 40.30959| 129| 242| 1231| 12| 1|
## | -7.473340| 43.78935| 140| 179| 931| 32| 1|
## | -6.575039| 42.91070| 84| 247| 1012| 39| 1|
## | -5.132756| 43.49572| 133| 190| 822| 40| 1|
## | -7.787378| 40.39362| 137| 247| 1143| 11| 1|
## | -4.941888| 43.35310| 128| 194| 843| 41| 1|
## | -7.621731| 40.34170| 101| 229| 1514| 15| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5940 -0.8312 -0.5389 0.8694 2.3942
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 23.2786567 4.2708072 5.451 5.02e-08 ***
## bio1 -0.0756324 0.0122819 -6.158 7.36e-10 ***
## bio7 -0.0480691 0.0085500 -5.622 1.89e-08 ***
## bio12 -0.0013621 0.0007445 -1.829 0.0673 .
## bio14 -0.0742556 0.0183522 -4.046 5.21e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 729.19 on 778 degrees of freedom
## Residual deviance: 634.11 on 774 degrees of freedom
## AIC: 959.93
##
## Number of Fisher Scoring iterations: 4
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 263
## n absences : 516
## AUC : 0.7402548
## cor : 0.3491264
## max TPR+TNR at : -0.05057615
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 263
## n absences : 10000
## AUC : 0.4897205
## cor : 0.01811699
## max TPR+TNR at : 0.3785007
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0.0005745427, 0.999671 (min, max)
##
##
##
## Notes:
##
## $separate.glms$cyreni
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9175fe23c8>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -5.321760| 40.35238| 85| 294| 625| 19| 1|
## | -3.821750| 40.87999| 94| 283| 558| 20| 1|
## | -4.014862| 40.73640| 91| 284| 555| 21| 1|
## | -4.107080| 40.71737| 91| 284| 555| 21| 1|
## | -4.179445| 40.64695| 101| 290| 487| 18| 1|
## | -4.127770| 40.67088| 91| 284| 555| 21| 1|
## | -3.939233| 40.80118| 94| 284| 568| 21| 1|
## | -4.120000| 40.60000| 119| 291| 452| 14| 1|
## | -5.180000| 40.32000| 89| 296| 646| 19| 1|
## | -5.770000| 40.39000| 110| 293| 574| 13| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7919 -0.5906 -0.4148 0.5079 1.9350
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -36.919128 15.723503 -2.348 0.0189 *
## bio1 -0.015420 0.039642 -0.389 0.6973
## bio7 0.105882 0.041494 2.552 0.0107 *
## bio12 0.008238 0.004293 1.919 0.0550 .
## bio14 0.220501 0.140888 1.565 0.1176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 210.72 on 205 degrees of freedom
## Residual deviance: 147.60 on 201 degrees of freedom
## AIC: 211.95
##
## Number of Fisher Scoring iterations: 4
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 76
## n absences : 130
## AUC : 0.8466599
## cor : 0.58633
## max TPR+TNR at : -0.04999943
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 76
## n absences : 10000
## AUC : 0.5956868
## cor : 0.03365248
## max TPR+TNR at : 0.6980342
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 2.2272e-09, 1 (min, max)
##
##
##
## Notes:
##
##
## $combined.glm
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9191adee18>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -5.171215| 43.06957| 78| 249| 917| 48| 1|
## | -6.036635| 43.02531| 76| 246| 1012| 46| 1|
## | -7.679727| 40.38852| 137| 247| 1143| 11| 1|
## | -7.790437| 40.30959| 129| 242| 1231| 12| 1|
## | -7.473340| 43.78935| 140| 179| 931| 32| 1|
## | -6.575039| 42.91070| 84| 247| 1012| 39| 1|
## | -5.132756| 43.49572| 133| 190| 822| 40| 1|
## | -7.787378| 40.39362| 137| 247| 1143| 11| 1|
## | -4.941888| 43.35310| 128| 194| 843| 41| 1|
## | -7.621731| 40.34170| 101| 229| 1514| 15| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5449 -0.8446 -0.5548 0.8645 2.3290
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 23.9228497 3.4964120 6.842 7.80e-12 ***
## bio1 -0.0793702 0.0097853 -8.111 5.01e-16 ***
## bio7 -0.0466212 0.0072213 -6.456 1.07e-10 ***
## bio12 -0.0015825 0.0006535 -2.421 0.0155 *
## bio14 -0.0839101 0.0154101 -5.445 5.18e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 939.91 on 984 degrees of freedom
## Residual deviance: 812.17 on 980 degrees of freedom
## AIC: 1197
##
## Number of Fisher Scoring iterations: 4
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 339
## n absences : 646
## AUC : 0.7493561
## cor : 0.3615897
## max TPR+TNR at : -0.06311968
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 339
## n absences : 10000
## AUC : 0.4916912
## cor : 0.02234313
## max TPR+TNR at : 0.3612526
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0.0004039292, 0.9996541 (min, max)
##
##
##
## Notes:
##
## $separate.aic
## [1] 1171.882
##
## $combined.aic
## [1] 1197.031araur <- combine.species(list(iberolacerta.clade$species$aranica, iberolacerta.clade$species$aurelioi))
# Says to keep cyreni seprate
ib.moses.4 <- moses.list(list(araur, iberolacerta.clade$species$bonnali), env = euro.worldclim)
ib.moses.4## $separate.glms
## $separate.glms$`aranica + aurelioi`
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9192e1d038>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | 1.572401| 42.99966| 105| 249| 871| 57| 1|
## | 1.103258| 42.77680| 41| 229| 1271| 86| 1|
## | 0.860000| 42.75000| 37| 228| 1294| 87| 1|
## | 0.860000| 42.84000| 82| 244| 1051| 72| 1|
## | -0.980000| 42.84000| 74| 243| 1096| 63| 1|
## | 0.980000| 42.75000| 37| 228| 1294| 87| 1|
## | 0.923200| 42.70100| 37| 228| 1294| 87| 1|
## | 0.804300| 42.60900| 33| 228| 1294| 86| 1|
## | 0.801100| 42.69900| 59| 237| 1166| 78| 1|
## | 0.798000| 42.78900| 59| 237| 1166| 78| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9830 -0.5561 -0.3186 -0.0681 3.1635
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -19.378453 15.964267 -1.214 0.225
## bio1 0.006335 0.031603 0.200 0.841
## bio7 0.043132 0.051251 0.842 0.400
## bio12 0.003369 0.004470 0.754 0.451
## bio14 0.074818 0.071899 1.041 0.298
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 116.449 on 212 degrees of freedom
## Residual deviance: 85.555 on 208 degrees of freedom
## AIC: 54.437
##
## Number of Fisher Scoring iterations: 6
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 42
## n absences : 171
## AUC : 0.8529657
## cor : 0.4181947
## max TPR+TNR at : 0.5931098
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 42
## n absences : 10000
## AUC : 0.7901738
## cor : 0.0703204
## max TPR+TNR at : 0.6440013
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 2.591377e-05, 0.9995366 (min, max)
##
##
##
## Notes:
##
## $separate.glms$bonnali
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9192d637e8>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | 0.663000| 42.68600| 51| 234| 1204| 80| 1|
## | 0.034285| 43.16230| 100| 245| 913| 56| 1|
## | 0.675000| 42.68600| 59| 237| 1166| 78| 1|
## | -0.445610| 42.98187| 51| 233| 1182| 74| 1|
## | 0.726000| 42.63300| 33| 228| 1294| 86| 1|
## | 0.981000| 42.66500| 25| 226| 1352| 90| 1|
## | 0.133631| 42.63613| 61| 238| 1053| 67| 1|
## | 0.860000| 42.66000| 25| 226| 1352| 90| 1|
## | 0.130000| 42.64000| 61| 238| 1053| 67| 1|
## | 0.990000| 42.48000| 59| 237| 1137| 76| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4979 -0.5890 -0.1963 0.5117 2.5442
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.53980 27.45079 -0.056 0.95527
## bio1 -0.10833 0.04191 -2.585 0.00974 **
## bio7 0.05973 0.09790 0.610 0.54180
## bio12 -0.01529 0.01334 -1.146 0.25162
## bio14 0.15187 0.13326 1.140 0.25442
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 155.26 on 172 degrees of freedom
## Residual deviance: 100.00 on 168 degrees of freedom
## AIC: 55.882
##
## Number of Fisher Scoring iterations: 6
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 56
## n absences : 117
## AUC : 0.8624847
## cor : 0.5327598
## max TPR+TNR at : 0.0129469
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 56
## n absences : 10000
## AUC : 0.6346929
## cor : 0.0475505
## max TPR+TNR at : 0.5031617
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 2.147658e-11, 0.9949861 (min, max)
##
##
##
## Notes:
##
##
## $combined.glm
##
##
## Formula: presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f91907e4d50>
##
##
## Data table (top ten lines):
##
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | 1.572401| 42.99966| 105| 249| 871| 57| 1|
## | 1.103258| 42.77680| 41| 229| 1271| 86| 1|
## | 0.860000| 42.75000| 37| 228| 1294| 87| 1|
## | 0.860000| 42.84000| 82| 244| 1051| 72| 1|
## | -0.980000| 42.84000| 74| 243| 1096| 63| 1|
## | 0.980000| 42.75000| 37| 228| 1294| 87| 1|
## | 0.923200| 42.70100| 37| 228| 1294| 87| 1|
## | 0.804300| 42.60900| 33| 228| 1294| 86| 1|
## | 0.801100| 42.69900| 59| 237| 1166| 78| 1|
## | 0.798000| 42.78900| 59| 237| 1166| 78| 1|
##
##
## Model:
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1,
## 2)], weights = weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2470 -0.5789 -0.3044 0.4633 3.5457
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -16.881063 13.144764 -1.284 0.1991
## bio1 -0.035117 0.019333 -1.816 0.0693 .
## bio7 0.059289 0.045317 1.308 0.1908
## bio12 0.002790 0.003887 0.718 0.4728
## bio14 0.030720 0.051901 0.592 0.5539
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 271.71 on 385 degrees of freedom
## Residual deviance: 189.43 on 381 degrees of freedom
## AIC: 103.24
##
## Number of Fisher Scoring iterations: 6
##
##
##
## Model fit (training data): class : ModelEvaluation
## n presences : 98
## n absences : 288
## AUC : 0.8552473
## cor : 0.4808458
## max TPR+TNR at : -0.1496999
##
##
## Environment space model fit (training data): class : ModelEvaluation
## n presences : 98
## n absences : 10000
## AUC : 0.8438204
## cor : 0.1341083
## max TPR+TNR at : 0.4625696
##
##
## Proportion of data wittheld for model fitting: FALSE
##
## Model fit (test data): [1] NA
##
##
## Environment space model fit (test data): [1] NA
##
##
## Suitability:
## class : RasterLayer
## dimensions : 54, 162, 8748 (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -10, 17, 39, 48 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 2.614036e-06, 0.9985111 (min, max)
##
##
##
## Notes:
##
## $separate.aic
## [1] 110.3187
##
## $combined.aic
## [1] 103.2383# Says to keep separate but not super convincingARC and AOC
Let’s start out by looking at range overlap as a function of time in Iberolacerta using raster maps:
ib.arc <- enmtools.aoc(iberolacerta.clade, overlap.source = "range", nreps = 10)
ib.arc##
##
## Age-Overlap Correlation test
##
## 10 replicates
##
## p values:
## (Intercept) empirical.df$age
## 0.1818182 0.1818182
And follow that up by looking at niche overlap using Bioclim models:
ib.aoc.bc <- enmtools.aoc(iberolacerta.clade, overlap.source = "bc", nreps = 10, env = euro.worldclim)